{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aade5bd-d68d-4485-9743-2da347bcb224",
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "Created on 2023-07-20\n",
    "势流理论圆柱绕流气动力分析 - 带环量修正\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import simps\n",
    "import pandas as pd\n",
    "\n",
    "# ======================== 物理模型参数类 ========================\n",
    "class FlowParameters:\n",
    "    \"\"\"定义流体力学计算参数\"\"\"\n",
    "    def __init__(self):\n",
    "        self.a = 1.0          # 圆柱半径 (m)\n",
    "        self.U = 1.0          # 来流速度 (m/s)\n",
    "        self.rho = 1.225      # 空气密度 (kg/m³)\n",
    "        self.Gamma_list = [    # 环量参数 (m²/s)\n",
    "            0, \n",
    "            2*np.pi*self.a*self.U, \n",
    "            4*np.pi*self.a*self.U\n",
    "        ]\n",
    "        self.theta = np.linspace(0, 2*np.pi, 1000)  # 角度离散化\n",
    "\n",
    "# ======================== 理论计算模块 ========================\n",
    "class PotentialFlowCalculator:\n",
    "    \"\"\"势流理论气动力计算器\"\"\"\n",
    "    @staticmethod\n",
    "    def pressure_coefficient(theta, Gamma=None, params=None):\n",
    "        \"\"\"\n",
    "        计算压力系数分布\n",
    "        参数：\n",
    "            theta : 方位角数组 (rad)\n",
    "            Gamma : 环量值 (m²/s)\n",
    "            params: FlowParameters实例\n",
    "        返回：\n",
    "            Cp : 压力系数数组\n",
    "        \"\"\"\n",
    "        a, U = params.a, params.U\n",
    "        \n",
    "        if Gamma is None or Gamma == 0:\n",
    "            # 无环量情况\n",
    "            return 1 - 4*np.sin(theta)**2\n",
    "        else:\n",
    "            # 含环量修正\n",
    "            term1 = -4*np.sin(theta)**2\n",
    "            term2 = -(2*Gamma/(np.pi*a*U)) * np.sin(theta)\n",
    "            term3 = -(Gamma/(2*np.pi*a*U))**2\n",
    "            return 1 + term1 + term2 + term3\n",
    "\n",
    "    @staticmethod\n",
    "    def theoretical_lift(Gamma, params):\n",
    "        \"\"\"库塔-茹科夫斯基升力理论值\"\"\"\n",
    "        return -Gamma / (params.a * params.U)\n",
    "\n",
    "# ======================== 数值积分模块 ========================\n",
    "class AerodynamicIntegrator:\n",
    "    \"\"\"气动力数值积分器\"\"\"\n",
    "    @staticmethod\n",
    "    def compute_forces(theta, Cp):\n",
    "        \"\"\"\n",
    "        计算无量纲气动力系数\n",
    "        参数：\n",
    "            theta : 方位角数组 (rad)\n",
    "            Cp    : 压力系数数组\n",
    "        返回：\n",
    "            CL : 升力系数\n",
    "            CD : 阻力系数\n",
    "        \"\"\"\n",
    "        integrand_lift = Cp * np.sin(theta)\n",
    "        integrand_drag = Cp * np.cos(theta)\n",
    "        \n",
    "        CL = (1/(2*np.pi)) * simps(integrand_lift, theta)\n",
    "        CD = (1/(2*np.pi)) * simps(integrand_drag, theta)\n",
    "        return np.round(CL, 6), np.round(CD, 6)  # 保留6位小数\n",
    "\n",
    "# ======================== 可视化模块 ========================\n",
    "class ResultVisualizer:\n",
    "    \"\"\"结果可视化工具\"\"\"\n",
    "    @staticmethod\n",
    "    def plot_pressure_distribution(theta, Cp_data, Gamma_list):\n",
    "        \"\"\"绘制压力系数分布曲线\"\"\"\n",
    "        plt.figure(figsize=(12, 6))\n",
    "        for i, (Gamma, Cp) in enumerate(Cp_data):\n",
    "            plt.subplot(1, len(Cp_data), i+1)\n",
    "            plt.plot(theta, Cp, 'r-', linewidth=2)\n",
    "            plt.title(f\"Γ = {Gamma:.1f} m²/s\", fontsize=12)\n",
    "            plt.xlabel(r'$\\theta$ (rad)', fontsize=10)\n",
    "            plt.ylabel('$C_p$', fontsize=10)\n",
    "            plt.grid(True, linestyle='--', alpha=0.7)\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "# ======================== 主程序 ========================\n",
    "if __name__ == \"__main__\":\n",
    "    # 初始化参数\n",
    "    params = FlowParameters()\n",
    "    calculator = PotentialFlowCalculator()\n",
    "    integrator = AerodynamicIntegrator()\n",
    "    \n",
    "    # 计算所有工况\n",
    "    results = []\n",
    "    Cp_data = []\n",
    "    for Gamma in params.Gamma_list:\n",
    "        # 计算压力系数\n",
    "        Cp = calculator.pressure_coefficient(params.theta, Gamma, params)\n",
    "        \n",
    "        # 计算气动力系数\n",
    "        CL, CD = integrator.compute_forces(params.theta, Cp)\n",
    "        theory_CL = calculator.theoretical_lift(Gamma, params)\n",
    "        \n",
    "        # 保存结果\n",
    "        results.append({\n",
    "            'Gamma': Gamma,\n",
    "            'Theory_CL': theory_CL,\n",
    "            'Numerical_CL': CL,\n",
    "            'Numerical_CD': CD,\n",
    "            'Relative_Error (%)': abs((CL - theory_CL)/theory_CL)*100 if theory_CL !=0 else 0\n",
    "        })\n",
    "        Cp_data.append( (Gamma, Cp) )\n",
    "    \n",
    "    # 结果表格输出\n",
    "    df = pd.DataFrame(results)\n",
    "    print(\"\\n\" + \"=\"*60)\n",
    "    print(\"气动力系数对比分析表\")\n",
    "    print(\"=\"*60)\n",
    "    print(df.round(4).to_markdown(index=False))\n",
    "    \n",
    "    # 可视化\n",
    "    ResultVisualizer.plot_pressure_distribution(params.theta, Cp_data, params.Gamma_list)\n",
    "\n",
    "    # 理论局限性分析\n",
    "    print(\"\\n\" + \"=\"*60)\n",
    "    print(\"理论模型局限性总结\")\n",
    "    print(\"=\"*60)\n",
    "    print(\"1. 无粘假设导致零阻力 (数值阻力: {:.1e}量级)\".format(df['Numerical_CD'].abs().max()))\n",
    "    print(\"2. 环量需人为指定，无法反映启动涡生成机制\")\n",
    "    print(\"3. 无法预测流动分离现象 (实际Re>40会出现卡门涡街)\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
